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The solvation of hydrophobic solutes in water is special because liquid and gas are 
almost at coexistence. In the common hypernetted chain approximation to integral 
equations, or equivalently in the homogenous reference fluid of molecular density 
functional theory, coexistence is not taken into account. Hydration structures and 
energies of nanometer-scale hydrophobic solutes are thus incorrect. In this article, we 
propose a bridge functional that corrects this thermodynamic inconsistency by intro¬ 
ducing a metastable gas phase for the homogeneous solvent. We show how this can 
be done by a third order expansion of the functional around the bulk liquid density 
that imposes the right pressure and the correct second order derivatives. Although 
this theory is not limited to water, we apply it to study hydrophobic solvation in 
water at room temperature and pressure and compare the results to all-atom simu¬ 
lations. The solvation free energy of small molecular solutes like n-alkanes and hard 
sphere solutes whose radii range from angstroms to nanometers is now in quantitative 
agreement with reference all atom simulations. The macroscopic liquid-gas surface 
tension predicted by the theory is comparable to experiments. This theory gives 
an alternative to the empirical hard sphere bridge correction used so far by several 
authors. 
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I. INTRODUCTION 


Implicit solvation techniques based on liquid-state theory such as integral equation theory 
in the interaction-sitd^ or molecular picturd^ or classical density functional theorj^^EI have 
proven to be successful for the computation of solvation properties. Those methods have 
shown to give thermodynamic and structural results that get closer and closer to all-atom 
simulations at a much lower numerical cost. A current challenge lies in the development 
and implementations of three-dimensional implicit solvation theories to describe molecular 
liquids and solutions. Recent developments in this direction have focused on Gaussian 
helcP theoretical approaches, or the 3D reference interaction site model (3D-RISM),*^*^ an 
appealing integral equation theory that has proven recently to be applicable to, e.g., structure 
prediction in complex biomolecular systems. Integral equations are, however, restricted 
by the choice of a closure relation, typically, Hypernetted Chain (HNC), Percus-Yevick or 
Kovalenko-Hirata. Despite their great potential, they remain difficult to control and improve, 
especially for arbitrary three-dimensional molecules, and they can prove difficult to converge. 


We have proposed recently a three dimensional formulation of molecular density func¬ 
tional theory (MDFT) in the homogeneous reference fluid approximation (HRF) to study 
solvatiorpl^. It has proven successful in studying solvation properties of solutes of arbitrary 
three-dimensional complexity embedded in various molecular solvents. However, when one 
comes to water, the HRF approximation fails even qualitatively to predict the solvation of 
large hydrophobic solute^. Such limitation can be explained by two essential features of 
water at ambient conditions that are not properly described by HRF functional. First, it is 
known that the solvation free energy of mesoscale apolar solutes can be modeled as the sum 
of a surface and volume termf^l. Por water, as well as for any solvent at room condition, the 
pressure is very low, inducing negligible PV term until large radiP^. Another key feature is 
that at ambient condition, water is close to its liquid-gas coexistence. As a consequence the 
solvation of big hydrophobic solutes may induce dewettin^®. 


In section [TT| we propose an extension of MDFT to introduce the liquid-gas coexistence 


of the solvent and to recover the correct pressure. Then, in section pITl we apply our theory 
to a model of water and compare the results of solvation of apolar solutes with reference 
all-atom Monte Carlo simulations (MC). 
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II. THEORY 


While the theory discussed here is generic to any classical density functional theory, it 
is described below in the framework of the MDFT for water introduced recentljA2*^. We 
start from the single point charge extended (SPC/E) model of watei^, that is, a model 
comprising one Lennard-Jones site and 3 partial charges. With MDFT, one computes the 
solvation free energy and the solvation structure of a solute of arbitrary shape that acts 
on the water density field through an external potential. This last quantity is the sum of 
an electrostatic vector field E{r) and a Lennard-Jones scalar field <hLj(^)^- In the general 
case, the functional of the density p{r, f2) depends upon the position r, and the molecular 
orientation of the (rigid) solvent molecule. There is no restriction on the molecular or 
chemical nature of the solvent molecule model, but to be rigid. In that particular model of 
water, p{r, f2) can be split into two distinct fields: the molecular density field n{r) coupled to 
‘I’Lj(n) and the polarization vector field P{r) coupled to E (r). These fields are themselves 
functionals of p{r, f2): 

n{r) = J p{r,n)dn, (1) 

P{r) = JJ p{r'— r')dr dfl. (2) 

where dfi denotes the integration over all molecular orientations. fi{r, H) = Jq <5(r— 

MSm(H))dM is the molecular polarization of a single water molecule at the origin of a carte¬ 
sian frame, gm and Sm are the charge and position of the m*^ solvent site. One should 
refers to Jeanmairet et aP for a complete description of MDFT for water. Without loss 
of generality, we will stick in what follows to solutes without partial charges, so that the 
polarization vector field is zero (P = 0 ). As a consequence, the free energy is, at dominant 
order, a functional of n(r) only. 

We now write the Helmholtz free energy functional, P[n\, that is the difference of the 
grand potential of the system containing the solute, 0, and without the solute, ©b. In this 
last case, the solvent is homogeneous at density n-Q (typically 1 g/cm^ for water): 

J-[n]=0[n]-0B. (3) 
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This leads to 


J^[n{r)] = keT 


n{r) In - — n{r) + n-Q 


Ub 


dr 


+ / n(r)<l>Lj(r)dr + j;^Jn(r)]. 


( 4 ) 


The terms of the right-hand side of Eqj^ corresponds to the usual decomposition^® into 
an ideal term accounting for information entropy, an external term accounting for the per¬ 
turbation by the solute through its external potential, and an excess term accounting for 
solvent-solvent correlations. This last, excess term, can be rewritten without additional 
approximation as 


J^e^cHr)] = 


knT 


An{r)c{r)An{r')drdr + 


_ ttHNC I 77 


b! 


(5) 


where r = ||r — r'||, An(r) = n(r) — hb, and c(r) is the direct correlation function of 
the homogeneous reference fluid at density riB- The first term thus corresponds to a series 
expansion in density of -Texc; around the density of the HRF, truncated at second order. 
Truncated information is put into an unknown bridge term, When = 0 , i.e. when 
we stick to the pure HRF approximation, Fqj^ can be shown to correspond to the HNC 
approximation of integral equation^. It will thus be called the HNC functional below. 
We suppose now that the correction can be expressed as a polynomial containing all terms 
of orders higher than 2 in An. Eqj^ can be used only if one knows the direct correlation 
function, c(r). In this article, we use an accurate direct correlation function of SPC/E water 
computed by Belloni et ah according to the methods discussed in refs. USEE! Fourier 
transform of the direct correlation function, c, is calculated by 

c(^ = —(6) 

1 -|- n^hooo {k) 

where h-ooo {k) is the first rotational invariant of the Fourier transform of the total correlation 
function h calculated by Puibassset and BellonP^. This function, as well as all higher 
rotational invariant components are obtained at short range by Monte Carlo sampling, and 
at higher range by integral equation closures, so that small k values are very accurate. 

The HNC functional has proved to be good enough for studying solvation in acetonitrile 
and in the Stockmayer fluid, but it exhibits wrong behaviors when coming to watei^^. To 
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improve the description of water we have proposed, as several other authord^^*^, an hard 
sphere bridge functional that consists in replacing all the unknown orders in An of the 
molecular fluid by the known ones of a hard sphere fluid of a diameter chosen on physical 
consideration^. Such a correction does improve the solvation of small molecular soluted^®^. 
However, the HNC functional or the functional with the hard sphere bridge, HNC+HSB, are 
not able to reproduce to date the solvation of hydrophobic solutes at both small and large 
length scales. This is an important discrepancy that originates from the fact that water at 
room conditions is close to liquid-gas coexistence and has a very low pressure, a fact that is 
impossible to account for consistently in HNCP^. 

We proposed recently a correction that imposes the essential physic^. It is based on the 
separation of the functional of Eq|^ with the hard sphere correction in a short range and a 
long range part. The long range part was then made compatible with the Van-der-Waals 
theory of phase coexistence at long range in a spirit similar to the hum-Chandler-Weeks 
theorj^. It introduces a coarse-grained density, similar in nature to weighted densities 
at the core of fundamental measure theories for hard sphere fluid^. We were then able to 
reproduce qualitatively the solvation of hydrophobic solutes at all length-scales. The surface 
tension was found too high, however, and the solvation structure of qualitative agreement 
only. It should be noted that the key role of the pressure of the fluid was not identified in this 
work: The pressure was consequently not explicitly considered as a control parameter even 
if this correction had an effect on the pressure. A functional that imposes the coexistence 
and the right pressure of the fluid is thus presented here. 

What we propose is an expression of Jq, that is cubic in An. There are two main motiva¬ 
tions to such an expression, (i) First, Evans et aP^and later Rickayzen and collaborator^^^^ 
showed that a series expansion of the functional at the quadratic order is thermodynami¬ 
cally inconsistent . In particular, the pressure of the homogeneous reference fluid predicted 
by the theory can be overestimated by orders of magnitude. For instance, for water, the 
HNC functional predicts a pressure of approximately 11450 bar instead of 1 bar. (u) Also, 
Rickaysen proposed to add the simplest cubic term to the series expansion of Wexc in den¬ 
sity, and showed it to be sufficient to overcome the thermodynamic inconsistency. Following 
Rickayzen’s prescriptions, a “three body” bridge term that is cubic in An, Wsb, is proposed. 
We give arguments on the form that should have Wsb for water, and how the addition of a 
physical constraint makes it a single parameter functional. 


5 






Instead of the simple three body expression of Rickayzen based on the overlap of hard 
bodies, we use here a rather different expression that is motivated by the fact that in water, 
tetrahedral order due to hydrogen bonding is lost in the HNC approximation and should be 
reinforced. Note that if the structuration discussed below is particular to water, the idea of 
including a three body term to improve the local structuration given by the HNC functional 
is relevant for any solvent. 

The three body functional should (i) enforce thermodynamic consistency, (ii) give back 
the local order due to A^-body interactions missed so far {N > 2), and (in) stay numerically 
efficient since it is our long-term goal to compete with other implicit methods like PCM 
(Polarizable Continuum Model)P^ that are much cruder but extremely useful. This last 
point may seem minor from a physical point of view; Nevertheless, to compute -Tsb, one 
should integrate over the whole instead of (with convolutions) for HNC: if it is not built 
efficiently, then it is useless. Consequently, in addition to physical motivation, the analytical 
form of the three body functional introduced here must allow efficient computation. 

We start from the idea of the coarse-grained model of tetracoordinated silicon by Stillinger 
and Webei®®, re-parameterized later for water by Molinero and Moord^. Their idea relies 
on an harmonic penalty to non-tetrahedral oxygen-oxygen-oxygen angles. In the MDFT 
framework, it leads to 


/3J'3B[n(r)] = 


X 


A 

2 


An(ri) 


An(r2)An(r3)/(ri2)/(ri3) 


/ ri2 • ^13 

V ri2ri^ 



dr2dr3 


dri. 


(7) 


with l3={k^T)~^. The dot product dehnes the cosine of the angle between three space points, 
and the quadratic term enforces a tetrahedral angle with 9q = 109.5°. The function / tunes 
the range of the three body interaction. As a source of local structuration of the fluid, it 
must be short-ranged and must vanish after few solvent radii, at distance r^ax- We propose 
as Molinero and Moore 


/(^) 


1 0 if r > Tj 


( 8 ) 


A is a dimensionless parameter modulating the strength of this oriented-bond term (hydrogen 
bond in case of water). The excess term in Eqj^is specihc to a given fluid and should thus 
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be parameterized once for all for the sake of consistency. We chose it so that one recovers 
the thermodynamic consistency and the correct pressnre of the bnlk liquid. 

The grand potential of a system of homogeneous fluid of volume V and pressure P is 
equal, by definition, to —PV. It is 0 in an empty system. Thus, one can deduce the pressure 
in the reference fluid by evaluating the functional of Eqj^ at zero densitj^ 


jr[n = 0] = e[n = 0] - 0B = PV, 

Using Eqj^for the functional without the three body term we get. 


n 


/^Phnc — fiB — 


(9) 


( 10 ) 


with c = 47r r^c{r)dr. With the three-body term: 


nt 


/^PsB =nB- + 9 


32?7,| 2 ^ 

71 A 


1 2 


/(r)r^dr 


L^o 


( 11 ) 


With Eqj^we find a pressure above 11450 bar for the HNC functional. Eqj^is used to fix 
the parameter A to have the desired pressure for the bulk fluid, i.e., 1 bar for water at room 
conditions. With this constraint, the three body functional has only one parameter left: the 
range of the interaction, r ma ^. Molinero and Moore determined a parameter r^ax = 4.3 A for 
their model. We kept the freedom of slightly varying around this value. An optimum 
value is found for 4.2 A. See below. 

The direct computation of the three-body function of Eqj^ cannot be performed because 
it requires a triple nested integration over the spacial coordinates. To accelerate the com¬ 
putation of this term we rewrite Eqj^ as: 


/?J3B[n(r)] = - 


where 


An(ri) I E UaisiriY + cos^ {9o) no{riY - 2 cos {9o) ni(ri) ■ ni(ri) j dri 

ya,l3£{x,y,z} 

( 12 ) 


n. 


*/3(n) = / f{ri2) 


«12/5i2 


An(r 2 )dr 2 , a,P e {x,y,z} 


12 


T12 


ni{ri) = / /(ri2)—An(r2)dr2, 

J ri2 

no{ri) = / /(ri2)An(r2)dr2. 


(13) 

(14) 

(15) 
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It can be seen that belongs to the general class of weighted functionals with one scalar 
weighted density, one vectorial one, and one second order, tensorial one. 

The derivation of the equivalence between Eqj^and Eqjl^as well as the hrst- and second- 
order functional derivatives that may be needed for minimizing Eqj^ are given in supple¬ 
mentary informat ioEp^. Convolution products of Eqs,]_3, ^and 1^ are evaluated efficiently 
in three dimensions using fast Fourier transforms (EFT). We typically use cubic boxes of 
35^ with space discretized by 5 grid nodes per A. Functional minimization of the total 
functional is typically reached within 15 to 20 iterations in a few tens of minutes on a single 
processor core at 2.4 GHz. 


III. RESULTS AND DISCUSSION 

Our goal is to predict the hydration structure and free energy of hydrophobic solutes from 
microscopic to macroscopic length scales. Hydration free energies of nanometric solutes are 
proportional to the surface of the solute. Since this behavior is due to the almost zero 
pressure of liquid water at room conditions, it is of prime importance to build a density 
functional that imposes the right pressure. First, we describe the parameterization of Eqj^ 
for capturing both the liquid-gas coexistence and the correct pressure. After that, we test 
the functional against hydration of various apolar solutes. 

To parametrize and test the three body term of Eqj^we hrst study small molecular apolar 
solutes. In Figj^ we show the solvation free energies of small alkane chains as computed by 
Monte Carlo simulations (MC)ESland by MDFT-HNC or MDFT-HNC+3B with = 4.3 A 
and 4.2 A. Within MDFT-HNC, the error in solvation free energy increases linearly with 
the size of the alkane, that is its number of carbons, shown here from methane to hexane. 
With the three-body excess functional, the relative error of MDFT with respect to Monte 
Carlo simulations is reduced by several orders. 

We hnd that r ma x = 4.2 A produces the optimal results, close to 4.3 A for Molinero 
and Moore and we stick to this value in the rest of the article. The remaining parameter 
A is chosen to impose the correct pressure in bulk water, P = 1 bar, from EqjTTj One 
hnds A = 38. We highlight that since the pressure of the fluid is now correct, the pressure 
correction term proposed by Sergiievskyi et aP^ is no longer required. 

The solvation free energy of n-alkanes into water is known to scale linearly with the 






Figure 1. Hydration free energies for the first six linear alkanes as calculated with MDFT-HNC 
and MDFT-HNC+3B, compared to Monte Carlo simulations by Ashbaugh et. al.l^. r^ax = 4.2 A 
(red squares) and 4.3 A (green circles) are shown for MDFT-HNC+3B. 



Fignre 2. Hydration free energy in SPC/E water for the hrst six linear alkanes as a function of the 
solvent accessible snrface area (SASA). Reference results from Monte Carlo are plotted as black 
circles. MDFT results with three-body corrections are in red squares. Linear regressions based on 
propane, butane, pentane and hexane are also plotted. 


molecular surface areaP^® 


AF = -k b, 


( 16 ) 


with A the solute area, 7 m the free energy per microscopic surface area and b an offset. 
Note that 7m is a microscopic equivalent to a surface tension, but is definitely different from 
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HNC 

HNC+3B 

MD 

Exp. 

7 m (J/(mol-A2)) 

340.59 

32.49 

30.53 

28.5 

b (kJ/mol) 

-15.02 

9.52 

7.81 

2.51 


Table I. Microscopic equivalent to the surface tension and offset from MC(^, from experiments^ 
and by MDFT-HNC and MDFT-HNC+3B. 


the macroscopic liquid-gas surface tension. Several definitions of A can be found in the 
litterature, that do not change any conclusion therein: We will use the solvent accessible 
surface area (SASA) of water in what follows, in order to be as comparable as possible 
with the results by Ashbaugh et al.l^SI. The evolution of the hydration free energy with 
respect to the solvent accessible surface area is plotted in Figj^ With the three-body excess 
functionnal, we now hnd the anticipated linear dependancy. The values and b given by 

are given in Table The value of 7 ^, is 


the linear regressions corresponding to equation 16 


in good agreement with both MC and experiments. We get an offset of approximately one 
keT with respect to MC. 


Now that parameters are fixed once for all, we show in Figj^the Helmholtz free energy 
of the homogenous systems as a function of the density at 300 K, as computed with the 
MDFT-HNC functional in dashed red and with the MDFT-HNC+3B functional in black. As 
discussed above, no second phase can appear in the system described with the MDFT-HNC 
functional since the free energy has only one minimum. Consequently, it can not capture 
liquid-gas coexistencd^. On the other hand, there are two minima of the Helmholtz free 
energy for the functional that includes the three-body bridge functional. A local minimum is 
found close to zero-density (“a gas phase”) with a free energy larger than the one of the global 
minimum corresponding to the density of the reference homogeneous fluid. The difference 
in Helmholtz free energy is of the order of 6 . 0 . 10 “^ kJ/A^, the homogeneous water we are 
describing is thus liquid and very close to liquid-gas coexistence. This physical feature is a 
keySSl to predict the solvation structure of large hydrophobic solutes of nanometer scale. 


To summarize: (i) the cost in free energy per unit volume for creating a cavity within 
the HNC (or HRF) formalism is several orders of magnitude too high, in relation to its 
overestimation of the pressure, (ii) the bridge functional that we propose corrects both the 
local order and the pressure, and it induces that the system is close to coexistence. The cost 
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Figure 3. Helmholtz free-energy of a homogeneous system of density n, see Eq. ne is the reference 
density one uses for the HNC functional. The insight is a focus on the first local minimum of the 
three-body corrected functional, HNC+3B. 

for creating a cavity within the MDFT-HNC+3B formalism is thus reduced to almost zero. 

We now focus on the solvation of hydrophobic solutes of atomic to nanoscale sizes. In 
their seminal works, Chandler and collaborator^^^^ studied by Monte Carlo simulations the 
hydration of hard spheres whose radii range from angstroms to nanometers. They observed a 
maximum in height of the first peak of the hard sphere - water radial distribution function at 
approximately 5 A. For radii larger than about 10 A, they also observed a slow convergence 
toward a plateau for the surface free energy. We compare the results by MDFT-HNC and 
MDFT-HNC+3B to those of Huang et al. in Figj^ Figj^and Figj^ One should keep in mind 
that MDFT results are approximatively 1000 times faster than explicit molecular dynamics 
or Monte Carlo simulations and that no other implicit solvent methods besides the LCW 
theory is able to reproduce these thermodynamic properties. 

In Figj^ we present the evolution of the height of the first peak of the hard sphere (HS) - 
water radial distribution function when the HS radius increases. This maximum corresponds 
to the most probable distance of molecules of the first solvation shell to the center of the 
hard sphere solute. The reference data by explicit methods are given in black?^. This height 
exhibits a peculiar maximum that is characteristic of the solvation of hydrophobic solutes in 
watei®. It tends toward unity for large radii. This behavior has been explained as follow: 
for small radii, the solvent can reorganize around the solute without losing solvent-solvent 
interactions, that is without losing too much cohesion: The increase of the height of the peak 


11 









n ^ ^ ^ ^ ^ ^ 

reference simulations 


Sy!»»>-yBraS Wl «... y « 



Distance to hard sphere solute (A) 


Figure 4. Radial distribution function of a hard sphere solute of growing radius in SPCE water at 
300 K from reference Monte Carlo simulation^ISl^ MDFT-HNC and MDFT-HNC+3B. 



Figure 5. Maxima of the radial distribution functions of hard sphere of different radii R. The MC 
simulations results of Huang et are the black circles, the ones obtained by MDFT-HNC+3B 
are the red crosses and the results of MDFT-HNC are the blue triangles. 

is due to an increase in packing of molecules at the surface of the sphere. For bigger radii, 
the perturbation is too high to keep the local structure unchanged: there is a loss of solvent- 
solvent interactions that has an energetic cost that limits the accumulation of molecules 
at the surface of the sphere and induces dewetting eventually. As a summary, when the 
perturbation stays small compared to solvent cohesion, the packing increases around the 
solute. Then, when the perturbation (the size of the solute) is unfavorable compared to 
solvent cohesion, solvent molecules stand back and the packing decreases. 
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MDFT-HNC fails to reproduce this behavior; as depicted in Figj^ there is no possible 
change of regime for the fluid. With the three body term, this change of regime can be found 
if the perturbation is able to make the fluid reach a state close to the second minimum. This 
is confirmed by Figj^ where MDFT-HNC+3B is in qualitative agreement with reference all 
atom simulations: The maximum of the radial distribution function is obtained around 
2.5 A, which is reasonable a value. The decrease is, however, too fast. 

In Figj^ we plot the solvation free energy of HS solutes per surface unit. We compile 
therein the results by MDFT-HNC, MDFT-HNC+3B and once again the reference all atom 
Monte Carlo simulations. MC shows a linear increase of the surface free energy for small 
radii, followed by a transition state, then followed by a plateau. This asymptotic value, 
reached at large HS radii corresponds to the surface tension of the fluid. At this regime, the 
solvation is thus driven by a sole surface term that corresponds at the microscopic level to 
the case where the loss of interaction between solvent particles is the prominent energetic 
term. MDFT-HNC is in agreement with the simulations only for very small radius (below 
2.5 A) but does not reproduce the plateau for bigger radius (> 10 A), this is consistent with 
the structural results, the transition between the two regimes is missed. Again, MDFT- 
HNC+3B is in good agreement with the simulations and the surface tension of SPC/E 
water estimated by Vega et ai^ is recovered. 

We can thus relate the decay of the maximum of the radial distribution function in Figj^ 
and the convergence to the plateau in Figj^ For structural and energetic properties, Monte 
Carlo simulations show a smooth transition between the two regimes described above, while 
MDFT-HNC+3B sharpens the transition: The three body term exacerbates the importance 
of the loss of attraction between solvent molecules. 

To conclude this section, (i) the structural properties obtained with MDFT-HNC+3B 
are improved with respect to MDFT-HNC since we recover the change in regime observed 
in MC at least qualitatively; (ii) this is also true for the the solvation free energy and 
this represents a considerable progress since MDFT-HNC predicts the wrong quantitative 
behavior, (in) The surface tension, that is related to the height of the saddle point in the 
free energy curve of Figj^ is correctly reproduced by MDFT-HNC+3B even though this is 
not explicitly controlled. 
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Figure 6. Solvation free energy for hard spheres of different radii R. The value of the liquid-vapor 
surface tension of SPC/E water at 300 K estimated by Vega et aP^ (63.6 mJ/m^) is shown in dotted 
black line. MDFT-HNC+3B gets the right behavior: two regimes at small then large HS radii, and 
the correct surface tension. 

IV. CONCLUSIONS 

In this paper, we propose to go beyond the usual quadratic expansion of the Gibbs free 
energy (or equivalently of the excess functional) around the homogeneous reference fluid 
within the molecular density functional theory framework. We thus go beyond the HNC 
approximation. MDFT-HNC+3B imposes a second local minimum to the Gibbs free energy 
of the system at low fluid density. The bridge functional that was proposed (i) enforces the 
tetrahedral order of water, (ii) recovers the close coexistence between gas and liquid states 
and their surface tension, and (in) is consistent with the experimental pressure of the fluid. 
It introduces one empirical parameter that we fix to parameterize over the solvation free 
energy of the first linear alkanes. It recovers the reference results of explicit simulations 
with a systematic offset of order keT. That is close to chemical accuracy, and is a clear 
improvement over MDFT-HNG. 

One advantage of this additional term with respect to previous work!^ is that (i) has 
a single empirical parameter, (ii) it does not require additional fields like coarse-grained 
densities, and (Hi) it makes the theory thermodynamically consistent. 

This bridge functional was used to study the solvation free energy of hard spheres whose 
radii range from angstroms to nanometers. Unlike MDFT-HNG, MDFT-HNG+3B recovers 
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the change of regime between a solvation governed by distortion of the solvent strnctnre 
and a solvation governed by a complete reorganization of the solvent. The free energy of 
solvation and the snrface tension are correct. Nevertheless, the transition stage is too sharp. 

Indeed, this points ont the necessity of fnrther improvements of the fnnctional. Other 
thermodynamic properties pertinent to hydrophobic solvation, snch as entropy, enthalpy, 
partial molar volnmes, temperatnre dependence, remain to be carefnlly tested too. 

Nnmerical efficiency is the very essence of implicit methods like MDFT. The bridge 
fnnctional introdnced therein wonld canse a dramatic increase of the nnmerical cost withont 
its rewriting in terms of fast Fonrier transforms. This is an important resnlt of this article. 
The nnmerical cost increase is at this stage of one order of magnitnde only with respect to 
MDFT-HNC. MDFT-HNC+3B is still two to three orders of magnitndes faster than explicit 
simnlations. 

Finally the solntes stndied here are all apolar and nentral, for the sake of clarity and 
pedagogy. The three-body fnnctional is bnilt to acconnt for short-range tetrahedral order in 
the solvent. It is similar in spirit to solnte-solvent corrections that were introdnced previonsly 
in the gronp to describe ions and H-bonded polar solnte^^^HSMl y^e think this will lead to 
a consistent fnnctional for water, valid for both hydrophobic and hydrophilic interactions. 
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